Take-Home Exercise 1

Published

January 30, 2023

Modified

February 12, 2023

1 Overview

1.1 Context

Water is an essential resource for people. The health of people depends on having access to clean water. It ensures peace and security, a clean environment, a sustainable economy, and a reduction in poverty. Yet more than 40% of the world’s population lacks access to enough clean water. According to UN-Water, 1.8 billion people will be residing in nations or areas with a severe water shortage by 2025. Water scarcity is a serious threat to several sectors, including food security. Around 70% of the world’s freshwater resources are used for agriculture.

1.2 Objectives

Discover the geographical distribution of functional and non-function water points and their co-locations if any in Osun State, Nigeria.

Exploratory Spatial Data Analysis (ESDA)

  • Derive kernel density maps of functional and non-functional water points. Using appropriate tmap functions

  • Display the kernel density maps on openstreetmap of Osub State, Nigeria.

  • Describe the spatial patterns revealed by the kernel density maps. Highlight the advantage of kernel density map over point map.

Second-order Spatial Point Patterns Analysis

With reference to the spatial point patterns observed in ESDA:

  • Formulate the null hypothesis and alternative hypothesis and select the confidence level.

  • Perform the test by using appropriate Second order spatial point patterns analysis technique.

  • With reference to the analysis results, draw statistical conclusions.

Spatial Correlation Analysis

  • Investigate if the spatial distribution of functional and non-functional water points are independent from each other.

  • Formulate the null hypothesis and alternative hypothesis and select the confidence level.

  • Perform the test by using appropriate Second order spatial point patterns analysis technique. With reference to the analysis results, draw statistical conclusions.

1.3 Data

Table 1: Data used
Type Name Format Description Source
Aspatial WPdx+ csv Locations of water points WPdx Global Data Repositories
Geospatial geoBoundaries SHP geoBoundaries data of Nigeria geoBoundaries

2 Importing and loading packages

For this exercise, we’ll be using the following packages:

  • sf: Manage and process vector-based geospatial data in R

  • spatstat: Perform 1st & 2nd spatial point patterns analysis + kernel density

  • raster: Reads, writes, manipulates, analyses, model of grid spatial data convert image output generate by spatstat into raster format

  • maptools: Convert Spatial objects into ppp format of spatstat

  • tmap: Plotting cartographic quality static point patterns maps or interactive maps by using leaflet API

  • tidyverse: a collection of functions, data, and documentation that extends the capabilities of base R.

  • funModeling: Plot charts for easier interpretations

We will be using p_load function of pacman package to install and load required packages.

pacman::p_load(maptools, sf, raster, spatstat, tmap, tidyverse, funModeling)

3 Spatial Data Wrangling

3.1 Data Geospatial

We will be using st_read function to read our geospatial data.

NGA <- st_read(dsn = "data/geospatial/", 
               layer = "nga_admbnda_adm2_osgof_20190417")%>%
  
st_transform(crs = 26392)
Reading layer `nga_admbnda_adm2_osgof_20190417' from data source 
  `C:\kt-x\is415-GAA\Take-Home_Ex\Take-Home_Ex01\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 774 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2.668534 ymin: 4.273007 xmax: 14.67882 ymax: 13.89442
Geodetic CRS:  WGS 84

Lets start by selecting specific fields that will be helpful in this exercise and select Osun state using the filter() function.

NGA <- NGA %>%
  select(c(3:4, 8:9)) %>%
  filter(ADM1_EN == "Osun")
glimpse(NGA)
Rows: 30
Columns: 5
$ ADM2_EN    <chr> "Aiyedade", "Aiyedire", "Atakumosa East", "Atakumosa West",…
$ ADM2_PCODE <chr> "NG030001", "NG030002", "NG030003", "NG030004", "NG030005",…
$ ADM1_EN    <chr> "Osun", "Osun", "Osun", "Osun", "Osun", "Osun", "Osun", "Os…
$ ADM1_PCODE <chr> "NG030", "NG030", "NG030", "NG030", "NG030", "NG030", "NG03…
$ geometry   <MULTIPOLYGON [m]> MULTIPOLYGON (((213526.6 34..., MULTIPOLYGON (…

3.2 Importing WPdx

We will be using read_csv() function to read the WPdx file and using filter function to only show Nigeria country data and only for Osun state.

wp_nga <- read_csv("data/aspatial/WPdx.csv") %>%
  filter(`#clean_country_name` == "Nigeria") %>%
  filter(`#clean_adm1` == "Osun")

3.2.1 Convert water point data into sf point features

wp_nga$Geometry = st_as_sfc(wp_nga$`New Georeferenced Column`)
wp_nga
# A tibble: 5,557 × 71
   row_id `#source`      #lat_…¹ #lon_…² #repo…³ #stat…⁴ #wate…⁵ #wate…⁶ #wate…⁷
    <dbl> <chr>            <dbl>   <dbl> <chr>   <chr>   <chr>   <chr>   <chr>  
 1 429123 GRID3             8.02    5.06 08/29/… Unknown <NA>    <NA>    Tapsta…
 2  70566 Federal Minis…    7.32    4.79 05/11/… No      Protec… Well    Mechan…
 3  70578 Federal Minis…    7.76    4.56 05/11/… No      Boreho… Well    Mechan…
 4  66401 Federal Minis…    8.03    4.64 04/30/… No      Boreho… Well    Mechan…
 5 422190 GRID3             7.87    4.88 08/29/… Unknown <NA>    <NA>    Tapsta…
 6 422064 GRID3             7.7     4.89 08/29/… Unknown <NA>    <NA>    Tapsta…
 7  65607 Federal Minis…    7.89    4.71 05/12/… No      Boreho… Well    Mechan…
 8  68989 Federal Minis…    7.51    4.27 05/07/… No      Boreho… Well    <NA>   
 9  67708 Federal Minis…    7.48    4.35 04/29/… Yes     Boreho… Well    Mechan…
10  66419 Federal Minis…    7.63    4.50 05/08/… Yes     Boreho… Well    Hand P…
# … with 5,547 more rows, 62 more variables: `#water_tech_category` <chr>,
#   `#facility_type` <chr>, `#clean_country_name` <chr>, `#clean_adm1` <chr>,
#   `#clean_adm2` <chr>, `#clean_adm3` <chr>, `#clean_adm4` <chr>,
#   `#install_year` <dbl>, `#installer` <chr>, `#rehab_year` <lgl>,
#   `#rehabilitator` <lgl>, `#management_clean` <chr>, `#status_clean` <chr>,
#   `#pay` <chr>, `#fecal_coliform_presence` <chr>,
#   `#fecal_coliform_value` <dbl>, `#subjective_quality` <chr>, …

3.2.2 Convert this df into sf object

wp_sf <- st_sf(wp_nga, crs=4326)
wp_sf
Simple feature collection with 5557 features and 70 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 4.032004 ymin: 7.060309 xmax: 5.06 ymax: 8.061898
Geodetic CRS:  WGS 84
# A tibble: 5,557 × 71
   row_id `#source`      #lat_…¹ #lon_…² #repo…³ #stat…⁴ #wate…⁵ #wate…⁶ #wate…⁷
 *  <dbl> <chr>            <dbl>   <dbl> <chr>   <chr>   <chr>   <chr>   <chr>  
 1 429123 GRID3             8.02    5.06 08/29/… Unknown <NA>    <NA>    Tapsta…
 2  70566 Federal Minis…    7.32    4.79 05/11/… No      Protec… Well    Mechan…
 3  70578 Federal Minis…    7.76    4.56 05/11/… No      Boreho… Well    Mechan…
 4  66401 Federal Minis…    8.03    4.64 04/30/… No      Boreho… Well    Mechan…
 5 422190 GRID3             7.87    4.88 08/29/… Unknown <NA>    <NA>    Tapsta…
 6 422064 GRID3             7.7     4.89 08/29/… Unknown <NA>    <NA>    Tapsta…
 7  65607 Federal Minis…    7.89    4.71 05/12/… No      Boreho… Well    Mechan…
 8  68989 Federal Minis…    7.51    4.27 05/07/… No      Boreho… Well    <NA>   
 9  67708 Federal Minis…    7.48    4.35 04/29/… Yes     Boreho… Well    Mechan…
10  66419 Federal Minis…    7.63    4.50 05/08/… Yes     Boreho… Well    Hand P…
# … with 5,547 more rows, 62 more variables: `#water_tech_category` <chr>,
#   `#facility_type` <chr>, `#clean_country_name` <chr>, `#clean_adm1` <chr>,
#   `#clean_adm2` <chr>, `#clean_adm3` <chr>, `#clean_adm4` <chr>,
#   `#install_year` <dbl>, `#installer` <chr>, `#rehab_year` <lgl>,
#   `#rehabilitator` <lgl>, `#management_clean` <chr>, `#status_clean` <chr>,
#   `#pay` <chr>, `#fecal_coliform_presence` <chr>,
#   `#fecal_coliform_value` <dbl>, `#subjective_quality` <chr>, …

3.2.3 Transforming into Nigeria projected coordinate system

wp_sf <- wp_sf %>%
  st_transform(crs = 26392)

3.2.4 EDA for Water Point Data

Next, we want to first have a visual of the water point data by creating a frequency chart. (Reference to In-class exercise 2)

freq(data = wp_sf,
     input = '#status_clean')

                     #status_clean frequency percentage cumulative_perc
1                       Functional      2319      41.73           41.73
2                   Non-Functional      2008      36.13           77.86
3                             <NA>       748      13.46           91.32
4      Functional but needs repair       248       4.46           95.78
5 Non-Functional due to dry season       151       2.72           98.50
6        Functional but not in use        63       1.13           99.63
7                        Abandoned        15       0.27           99.90
8         Abandoned/Decommissioned         5       0.09          100.00

Seems like there are several different categories of water points. With reference to Objectives, we will be focusing on Functional and Non-Functional while still taking into account of the unknown water points. But first, let’s rename this spatial frame.

wp_sf_nga <- wp_sf %>% 
  rename(status_clean = '#status_clean') %>%
  select(status_clean) %>%
  mutate(status_clean = replace_na(
    status_clean, "unknown"))

After that, we filter and assign the wp_functional

wp_functional <- wp_sf_nga %>%
  filter(status_clean %in%
           c("Functional",
             "Functional but not in use",
             "Functional but needs repair"))
freq(data = wp_functional,
     input = 'status_clean')

                 status_clean frequency percentage cumulative_perc
1                  Functional      2319      88.17           88.17
2 Functional but needs repair       248       9.43           97.60
3   Functional but not in use        63       2.40          100.00

Sames goes for the Non-Functional

wp_nonfunctional <- wp_sf_nga %>%
  filter(status_clean %in%
           c("Abandoned/Decommissioned",
             "Abandoned",
             "Non-Functional due to dry season",
             "Non-Functional",
             "Non functional due to dry season"))
freq(data = wp_nonfunctional,
     input = 'status_clean')

                      status_clean frequency percentage cumulative_perc
1                   Non-Functional      2008      92.15           92.15
2 Non-Functional due to dry season       151       6.93           99.08
3                        Abandoned        15       0.69           99.77
4         Abandoned/Decommissioned         5       0.23          100.00

Sames goes for the Unknown

wp_unknown <- wp_sf_nga %>%
  filter(status_clean == "unknown")

freq(data = wp_unknown,
     input = 'status_clean')

  status_clean frequency percentage cumulative_perc
1      unknown       748        100             100

3.2.5 Performing Point-in-Polygon Count

Next, we want to find out the number of total, functional, nonfunctional and unknown water points in each LGA. First, it identifies the functional water points in each LGA by using st_intersects() of sf package. Next, length() is used to calculate the number of functional water points that fall inside each LGA.

NGA_wp <- NGA %>% 
  mutate(`total_wp` = lengths(
    st_intersects(NGA, wp_sf_nga))) %>%
  mutate(`wp_functional` = lengths(
    st_intersects(NGA, wp_functional))) %>%
  mutate(`wp_nonfunctional` = lengths(
    st_intersects(NGA, wp_nonfunctional))) %>%
  mutate(`wp_unknown` = lengths(
    st_intersects(NGA, wp_unknown)))

4 Geospatial Data wrangling

Moving to data wrangling on geospatial data, first, we need to convert simple feature data frame to sp spatial class.

4.1 Converting sf data frames to sp’s Spatial* class

Convert to sp object/class using as_Spatial(). Take a look at the following 2 tabs (wp_sc, NGA_sc), notice the properties is “SpatialPolygonsDataFrame”.

wp_sc_functional <- as_Spatial(wp_functional)
wp_sc_nonfunctional <- as_Spatial(wp_nonfunctional)
NGA_sc <- as_Spatial(NGA)
wp_sc_functional
class       : SpatialPointsDataFrame 
features    : 2630 
extent      : 177285.9, 290751, 343128.1, 450859.7  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 
variables   : 1
names       :              status_clean 
min values  :                Functional 
max values  : Functional but not in use 
wp_sc_nonfunctional
class       : SpatialPointsDataFrame 
features    : 2179 
extent      : 180539, 290616, 340054.1, 450780.1  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 
variables   : 1
names       :                     status_clean 
min values  :                        Abandoned 
max values  : Non-Functional due to dry season 
NGA_sc
class       : SpatialPolygonsDataFrame 
features    : 30 
extent      : 176503.2, 291043.8, 331434.7, 454520.1  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 
variables   : 4
names       :  ADM2_EN, ADM2_PCODE, ADM1_EN, ADM1_PCODE 
min values  : Aiyedade,   NG030001,    Osun,      NG030 
max values  :   Osogbo,   NG030030,    Osun,      NG030 

4.2 Converting the Spatial* class into generic sp format

WHY? Because spatstat requires the analytical data in ppp object form. There is no direct way to convert a Spatial* classes into ppp object. We need to convert the Spatial classes* into Spatial object first.

Starting with Converting the Spatial* classes into generic sp objects. Take a look at the following 2 tabs (wp_sc, NGA_sc), notice the properties is “SpatialPoints” and “SpatialPolygons” respectively.

#wp_sp <- as(wp_sc, "SpatialPoints")
wp_sp_functional <- as(wp_sc_functional, "SpatialPoints")
wp_sp_nonfunctional <- as(wp_sc_nonfunctional, "SpatialPoints")
nga_sp <- as(NGA_sc, "SpatialPolygons")
#wp_sp
wp_sp_functional
class       : SpatialPoints 
features    : 2630 
extent      : 177285.9, 290751, 343128.1, 450859.7  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 
wp_sp_nonfunctional
class       : SpatialPoints 
features    : 2179 
extent      : 180539, 290616, 340054.1, 450780.1  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 
nga_sp
class       : SpatialPolygons 
features    : 30 
extent      : 176503.2, 291043.8, 331434.7, 454520.1  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=4 +lon_0=8.5 +k=0.99975 +x_0=670553.98 +y_0=0 +a=6378249.145 +rf=293.465 +towgs84=-92,-93,122,0,0,0,0 +units=m +no_defs 

4.3 Converting the generic sp format into spatstat’s ppp format

Next, we will use as.ppp() function of spatstat to convert the spatial data into spatstat’s ppp object format.

#wp_ppp <- as(wp_sp, "ppp")
wp_ppp_functional <- as(wp_sp_functional, "ppp")
wp_ppp_nonfunctional <- as(wp_sp_nonfunctional, "ppp")

#wp_ppp
wp_ppp_functional
Planar point pattern: 2630 points
window: rectangle = [177285.9, 290750.96] x [343128.1, 450859.7] units
wp_ppp_nonfunctional
Planar point pattern: 2179 points
window: rectangle = [180538.96, 290616] x [340054.1, 450780.1] units
#plot(wp_ppp)

par(mfrow=c(1,2))
plot(wp_ppp_functional)
plot(wp_ppp_nonfunctional)

4.3.1 Check for duplicates

#any(duplicated(wp_ppp))

any(duplicated(wp_ppp_functional))
[1] FALSE
any(duplicated(wp_ppp_nonfunctional))
[1] FALSE

If so, check for the number of duplicates

sum(multiplicity(wp_ppp_functional) > 1)
[1] 0
sum(multiplicity(wp_ppp_nonfunctional) > 1)
[1] 0

To resolve this problem, we will be using the jittering approach, which will add a small perturbation to the duplicate points so that they do not occupy the exact same space.

#wp_ppp_jit <- rjitter(wp_ppp, 
#                      retry=TRUE,
#                      nsim=1, 
#                      drop=TRUE)

wp_ppp_functional_jit <- rjitter(wp_ppp_functional, 
                      retry=TRUE,
                      nsim=1, 
                      drop=TRUE)
wp_ppp_nonfunctional_jit <- rjitter(wp_ppp_nonfunctional, 
                      retry=TRUE,
                      nsim=1, 
                      drop=TRUE)

Then, check again for duplicates

any(duplicated(wp_ppp_functional_jit))
[1] FALSE
any(duplicated(wp_ppp_nonfunctional_jit))
[1] FALSE

4.4 Create owin object

owin object is designed to represent this polygonal region. We will be using to convert nigeria Spatial Polygon object into owin object of spatstat.

nga_owin <- as(nga_sp, "owin")
plot(nga_owin)

summary(nga_owin)
Window: polygonal boundary
30 separate polygons (no holes)
            vertices      area relative.area
polygon 1        204 766084000       0.08870
polygon 2         81 304399000       0.03520
polygon 3         97 465688000       0.05390
polygon 4        124 373051000       0.04320
polygon 5         60 149473000       0.01730
polygon 6         84 144820000       0.01680
polygon 7         50 102243000       0.01180
polygon 8         72 216002000       0.02500
polygon 9        112 269897000       0.03130
polygon 10       125 365142000       0.04230
polygon 11        83 111191000       0.01290
polygon 12       126 192557000       0.02230
polygon 13       219 904397000       0.10500
polygon 14       174 741131000       0.08580
polygon 15        81 138742000       0.01610
polygon 16        65 119452000       0.01380
polygon 17        90 280205000       0.03240
polygon 18        69  69814600       0.00808
polygon 19        69  42727500       0.00495
polygon 20        49  30458800       0.00353
polygon 21        62 263505000       0.03050
polygon 22        93 438930000       0.05080
polygon 23        87 274127000       0.03170
polygon 24       105 509979000       0.05910
polygon 25        98 292058000       0.03380
polygon 26        64 327765000       0.03800
polygon 27       133 108945000       0.01260
polygon 28       122 462169000       0.05350
polygon 29        94 109715000       0.01270
polygon 30        95  61239800       0.00709
enclosing rectangle: [176503.22, 291043.82] x [331434.7, 454520.1] units
                     (114500 x 123100 units)
Window area = 8635910000 square units
Fraction of frame area: 0.613

4.4.1 Combining point events object and owin object

#wpNGA_ppp = wp_ppp[nga_owin]

wpNGA_ppp_functional = wp_ppp_functional[nga_owin]
wpNGA_ppp_nonfunctional= wp_ppp_nonfunctional[nga_owin]
#plot(wpNGA_ppp)

par(mfrow=c(1,2))
plot(wpNGA_ppp_functional)
plot(wpNGA_ppp_nonfunctional)

4.4.2 tmap plots

We can further plot our water points (functional & non-functional) using tmap().

Lets put both functional and non functional water points together. Also, Set the base map to be “OpenStreetMap”.

tmap_mode("view")
tm_basemap("OpenStreetMap") +
tm_shape(NGA_wp) + 
  tm_polygons() + 
tm_shape(wp_functional) + 
  tm_dots(col = "status_clean", 
          size = 0.01, 
          border.col = "black", 
          border.lwd = 0.5) +
tm_shape(wp_nonfunctional) + 
  tm_dots(col = "status_clean", 
          size = 0.01, 
          border.col = "black", 
          border.lwd = 0.5) +
  tm_view(set.zoom.limits = c(8, 16))

From this point map, even with the different colours points, it is difficult to intrepet the clusters of the functional and non-functional water points. Lets go to the next tab to look at only functional points.

tmap_mode("view")
tm_basemap("OpenStreetMap") +
tm_shape(NGA_wp) + 
  tm_polygons() + 
tm_shape(wp_functional) + 
  tm_dots(col = "status_clean", 
          size = 0.01, 
          border.col = "black", 
          border.lwd = 0.5) +
  tm_view(set.zoom.limits = c(8, 16))

At a macro view (zoom out), there are A LOT of functional water points. While zooming in, it looks like there’s some clusters of water points such as the border between Ife Central and Ife East. Lets take a look at the next tab, wp_nonfunctional.

tmap_mode("view")
tm_basemap("OpenStreetMap") +
tm_shape(NGA_wp) + 
  tm_polygons() + 
tm_shape(wp_nonfunctional) + 
  tm_dots(col = "status_clean", 
          size = 0.01, 
          border.col = "black", 
          border.lwd = 0.5) +
  tm_view(set.zoom.limits = c(8, 16))

Again, at a macro view (zoom out), there are A LOT of functional water points. While zooming in, it looks like there’s some clusters of water points such as the border between Ife Central and Ife East, which seems like the same cluster as shown in wp_functional. This highlights that point map might not suitable in identifying the clusters due to the overlapping points and categorisation. Thus, lets proceed onto the First-Order of Spatial Point Patterns Analysis to plot Kernel density maps which could provide more meaningful insights.


5 First-order Spatial Point Patterns Analysis

Now that we have the point maps from the previous section, moving forward, lets see how would the kernel density maps would display the water points in Osun.

5.1 Kernel Density Estimation

In this section, We will be computing Functional and Non-Functional water points in Osun.

5.1.1 Computing kernel density estimation using automatic bandwidth selection method

We will be using the bw.diggle() method than bw.ppl() because it is more suitable for this exercise to detect single tight cluster in the midst of random noise which was observed earlier section in the point maps, theres seem to be clusters, for example, along the border between Ife Central and Ife East.

kde_wpNGA_functional.bw <- density(wpNGA_ppp_functional,
                          sigma=bw.diggle,
                          edge=TRUE,
                          kernel="gaussian") 

kde_wpNGA_nonfunctional.bw <- density(wpNGA_ppp_nonfunctional,
                          sigma=bw.diggle,
                          edge=TRUE,
                          kernel="gaussian") 

par(mfrow=c(1,2))
plot(kde_wpNGA_functional.bw)
plot(kde_wpNGA_nonfunctional.bw)

It looks the unit of measurement is of the default value in meter. That’s why the density values computed is in “number of points per square meter”.

Before moving to the next section, let’s retrieve the bandwidth used to compute the KDE layer.

#bw <- bw.diggle(wpNGA_ppp)
#bw
bw_functional <- bw.diggle(wpNGA_ppp_functional)
bw_functional
   sigma 
252.1687 
bw_nonfunctional <- bw.diggle(wpNGA_ppp_nonfunctional)
bw_nonfunctional
   sigma 
308.2061 

5.1.2 Rescalling KDE values

wpNGA_ppp_functional.km <- rescale(wpNGA_ppp_functional, 1000, "km")
wpNGA_ppp_nonfunctional.km <- rescale(wpNGA_ppp_nonfunctional, 1000, "km")
#kde_wpNGA.bw <- density(wpNGA_ppp.km, sigma=bw.diggle, edge=TRUE, #kernel="gaussian")
#plot(kde_wpNGA.bw)


kde_wpNGA_functional.bw <- density(wpNGA_ppp_functional.km,
                          sigma=bw.diggle,
                          edge=TRUE,
                          kernel="gaussian") 

kde_wpNGA_nonfunctional.bw <- density(wpNGA_ppp_nonfunctional.km,
                          sigma=bw.diggle,
                          edge=TRUE,
                          kernel="gaussian") 

par(mfrow=c(1,2))
plot(kde_wpNGA_functional.bw)
plot(kde_wpNGA_nonfunctional.bw)

Unlike the point map, from the above KDE maps, it highlight the signifcant clusters and much easier to identify the clusters of functional and non-functional water points clusters.

5.2 Adaptive KDE

5.2.1 Computing KDE by using adaptive bandwidth

Fixed bandwidth method is very sensitive to highly skew distribution of spatial point patterns over geographical units for example urban versus rural. Hence, for this exercise, we will be using adaptive bandwidth instead.

Lets look at functional wp in adaptive bandwidth then move on the next tab for the non-functional.

#kde_wpNGA_adaptive <- adaptive.density(wpNGA_ppp.km, method="kernel")
#plot(kde_wpNGA_adaptive)

par(mfrow=c(1,2))

kde_wpNGA_functional_adaptive <- adaptive.density(wpNGA_ppp_functional.km, method="kernel")
plot(kde_wpNGA_functional_adaptive)
plot(kde_wpNGA_functional.bw)

kde_wpNGA_nonfunctional_adaptive <- adaptive.density(wpNGA_ppp_nonfunctional.km, method="kernel")
plot(kde_wpNGA_nonfunctional_adaptive)

plot(kde_wpNGA_nonfunctional.bw)

OBSERVATION Well it seems like theres isnt signficant changes to the functional BUT if u take a closer look at the non-functional map, some points are now lower in density with this adaptive bandwidth.

5.2.2 Converting KDE output into grid object

To make the KDE output is suitable for mapping purposes.

#gridded_kde_wpNGA_bw <- as.SpatialGridDataFrame.im(kde_wpNGA.bw)
#spplot(gridded_kde_wpNGA_bw)

gridded_kde_wpNGA_func_bw <- as.SpatialGridDataFrame.im(kde_wpNGA_functional.bw)
spplot(gridded_kde_wpNGA_func_bw)

gridded_kde_wpNGA_nfunc_bw <- as.SpatialGridDataFrame.im(kde_wpNGA_nonfunctional.bw)
spplot(gridded_kde_wpNGA_nfunc_bw)

5.2.2.1 Converting gridded output into raster

Next, we will convert the gridded kernal density objects into RasterLayer object by using raster() of raster package.

#kde_wpNGA_bw_raster <- raster(gridded_kde_wpNGA_bw)
#kde_wpNGA_bw_raster

kde_wpNGA_func_bw_raster <- raster(gridded_kde_wpNGA_func_bw)
kde_wpNGA_func_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.8948485, 0.9616045  (x, y)
extent     : 176.5032, 291.0438, 331.4347, 454.5201  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : v 
values     : -5.092436e-15, 25.49435  (min, max)
kde_wpNGA_nfunc_bw_raster <- raster(gridded_kde_wpNGA_nfunc_bw)
kde_wpNGA_nfunc_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.8948485, 0.9616045  (x, y)
extent     : 176.5032, 291.0438, 331.4347, 454.5201  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : v 
values     : -3.925434e-15, 20.49412  (min, max)

WHY? - crs property is NA during the convert

#projection(kde_wpNGA_bw_raster) <- CRS("+init=EPSG:3414")
#kde_wpNGA_bw_raster

projection(kde_wpNGA_func_bw_raster) <- CRS("+init=EPSG:3414")
kde_wpNGA_func_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.8948485, 0.9616045  (x, y)
extent     : 176.5032, 291.0438, 331.4347, 454.5201  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : v 
values     : -5.092436e-15, 25.49435  (min, max)
projection(kde_wpNGA_nfunc_bw_raster) <- CRS("+init=EPSG:3414")
kde_wpNGA_nfunc_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.8948485, 0.9616045  (x, y)
extent     : 176.5032, 291.0438, 331.4347, 454.5201  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : v 
values     : -3.925434e-15, 20.49412  (min, max)

Lets display the raster in cartographic quality map using tmap package.

tm_shape(kde_wpNGA_func_bw_raster) + 
  tm_raster("v") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE)
tm_shape(kde_wpNGA_nfunc_bw_raster) + 
  tm_raster("v") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE)

5.2.3 Spatial Point Patterns using KDE

5.2.3.1 Extract Osun

In this section, we will be comparing KDE of water points at Osun as the study area.

#osun = NGA_sc[wp_sc$clean_adm1 == "Osun",]
osun = NGA_sc

plot(osun, main = "Osun")

5.2.3.2 Converting the spatial point data frame into generic sp format

osun_sp = as(osun, "SpatialPolygons")

5.2.3.3 Creating owin object

osun_owin = as(osun_sp, "owin")

5.2.3.4 Combining wp points and Osun (Study Area)

wp_osun_ppp_func = wp_ppp_functional_jit[osun_owin]
wp_osun_ppp_nfunc = wp_ppp_nonfunctional_jit[osun_owin]

Next, rescale() function is used to trasnform the unit of measurement from metre to kilometre.

wp_osun_ppp_func.km = rescale(wp_osun_ppp_func, 1000, "km")
wp_osun_ppp_nfunc.km = rescale(wp_osun_ppp_nfunc, 1000, "km")

After that, we do the plotting

par(mfrow=c(1,2))
plot(wp_osun_ppp_func.km, main="Osun Functional")
plot(wp_osun_ppp_nfunc.km, main="Osun Non-Functional")

5.2.3.5 Computing KDE

The code chunk below will be used to compute the KDE of these four planning area. bw.diggle method is used to derive the bandwidth of each

par(mfrow=c(1,2))
plot(density(wp_osun_ppp_func.km,
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Osun Functional")

plot(density(wp_osun_ppp_nfunc.km,
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Osun Non-Functional")

6 Second-order Spatial Point Patterns

6.1 Analysing Spatial Point Process Using G-Function

In this section, the G function measures the distances between any two events and their respective nearest neighbours. Using spatstat package’s Gest() and envelope() to run a Monte Carlo simulation test, which is used to predict the probability of a range of outcomes when the possibility of random variables is present.

6.1.1 1. Computing G-function estimation (Functional)

Lets start with the functional.

G_func_osun = Gest(wp_osun_ppp_func, correction = "border")
plot(G_func_osun, xlim=c(0,500))

6.1.2 2. Performing Complete Spatial Randomness Test (Functional)

Next, we will be conducting a hypothesis test to confirm the observed spatial patterns above. The hypothesis and test are as follows:

Ho = The distribution of non-functional water points in Osun, Nigeria are randomly distributed.

H1= The distribution of non-functional water points in Osun, Nigeria are NOT randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.05.

Monte Carlo test with G-function

G_func_osun.csr <- envelope(wp_osun_ppp_func, Gest, nsim = 39)
Generating 39 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,  39.

Done.
plot(G_func_osun.csr)

Conclusion From func_osun.csr, it is observed that at 95% confidence interval, the G(r) far above the G(theo) and the envelope. This implies that the functional water points in Osun, Nigeria are clustered. Thus, we can reject the null hypothesis that the distribution of functional water points are randomly distributed.

6.1.3 3. Computing G-function estimation (Non-Functional)

Next, Non-functional

G_nfunc_osun = Gest(wp_osun_ppp_nfunc, correction = "border")
plot(G_nfunc_osun, xlim=c(0,500))

6.1.4 4. Performing Complete Spatial Randomness Test (Non-Functional)

Next, we will be conducting a hypothesis test to confirm the observed spatial patterns above. The hypothesis and test are as follows:

Ho = The distribution of non-functional water points in Osun, Nigeria are randomly distributed.

H1= The distribution of non-functional water points in Osun, Nigeria are NOT randomly distributed.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.05.

Monte Carlo test with G-function

G_nfunc_osun.csr <- envelope(wp_osun_ppp_nfunc, Gest, nsim = 39)
Generating 39 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,  39.

Done.
plot(G_nfunc_osun.csr)

Conclusion From nfunc_osun.csr, it is observed that at 95% confidence interval, the G(r) far above the G(theo) and the envelope. This implies that the non-functional water points in Osun, Nigeria are clustered. Thus, we can reject the null hypothesis that the distribution of non-functional water points are randomly distributed.

7 Analysing Spatial Point Process Using L-Function

Lets also investigate whether the spatial distribution of functional and non-functional water points are independent from each other.

7.0.1 Computing G-function estimation

# L_func_osun.csr = Lest(wp_osun_ppp_func, correction = "Ripley")
# plot(L_func_osun.csr, . -r ~ r, 
#      ylab= "L(d)-r", xlab = "d(m)")

7.0.2 Performing Complete Spatial Randomness Test

Again, we will be conducting a hypothesis test to confirm the observed spatial patterns above. The hypothesis and test are as follows:

Ho = The distribution of functional water points in Osun, Nigeria are spatially independent.

H1= The distribution of functional water points in Osun, Nigeria are NOT spatially independent.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.05.

Monte Carlo test with L-function

#L_func_osun.csr <- envelope(wp_osun_ppp_func, Lest, nsim = 39, rank = 1, glocal=TRUE)

7.0.3 Plot the result from L-function

#plot(L_func_osun.csr, . - r ~ r, xlab="d", ylab="L(d)-r")

CONCLUSION From the plot above, the empirical k-cross line is far outside from the envelope of the 95% confidence level, and for that scale, we reject the null hypothesis. Moreover, it is not even within the envelope in all distance(m). We can infer that the distribution of functional water points in Osun, Nigeria are NOT spatially independent.

  • NOTE: erm I was unable to run the output in time for this section for submission..so i had to comment these codes out in order to render

7.0.4 Computing G-function estimation

# L_nfunc_osun.csr = Lest(wp_osun_ppp_nfunc, correction = "Ripley")
# plot(L_nfunc_osun.csr, . -r ~ r, ylab= "L(d)-r", xlab = "d(m)")

7.0.5 Performing Complete Spatial Randomness Test

Again, we will be conducting a hypothesis test to confirm the observed spatial patterns above. The hypothesis and test are as follows:

Ho = The distribution of non-functional water points in Osun, Nigeria are spatially independent.

H1= The distribution of non-functional water points in Osun, Nigeria are spatially independent.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.05.

Monte Carlo test with L-function

#L_nfunc_osun.csr <- envelope(wp_osun_ppp_nfunc, Lest, nsim = 39, rank = 1, glocal=TRUE)

7.0.6 Plot the result from L-function

#plot(L_nfunc_osun.csr, . - r ~ r, xlab="d", ylab="L(d)-r")

CONCLUSION

8 Acknowledgement

Thank you Prof Kam for our IS415 Geospatial Analytics and Applications course materials & resources